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LONG-TERM  GOALS 

Our  long-term  goal  is  to  implement  an  operational  high-resolution  atmospheric  data  assimilation  and 
prediction  system  and  to  use  it  for  daily  weather  forecasting.  To  date  we  have  worked  on  several 
operational  and  scientific  aspects  of  the  problem  using  the  Coupled  Ocean/ Atmosphere  Mesoscale 
Prediction  System  (COAMPS®  ^);  (a)  setting  up  stable  large  scale  data  transfer  capabilities  to  facilitate 
COAMPS  runs  2-4  times  per  day  in  an  operational  manner  for  different  geographical  regions  and  using 
different  nested  grid  configurations;  (b)  perform  research  on  the  MPI  scalability  of  COAMPS  on 
selected  computer  architectures  and  to  optimize  the  code  to  take  advantage  of  the  vector  capabilities  of 
the  Cray  XI  and  the  massively  parallel  features  of  our  Linux  cluster;  (c)  identify  and  understand  the 
uncertainties  in  high-resolution  NWP  forecast  and  their  impact  on  severe  weather,  such  as  extreme 
rainfall,  and  to  develop  model  metrics  appropriate  to  mesoscale  weather  phenomena;  and  (d)  improve 
our  knowledge  of  observational  error  characteristics  for  spatially  correlated  data  and  develop  the 
numerical  schemes  capable  of  assimilating  these  types  of  observations. 


COAMPS®  is  a  registered  trademark  of  the  Naval  Research  Laboratory. 
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OBJECTIVES 


The  objectives  of  this  project  are  to:  (a)  implement  an  operational  data  feed  from  the  Navy  Operational 
Global  Atmospheric  Prediction  System  (NOGAPS),  and  implement  a  semi-operational  version  of 
COAMPS  at  the  Interdisciplinary  Centre  for  Mathematical  and  Computational  Modelling  (ICM), 
Warsaw  University;  (b)  validate  COAMPS  model  performance  through  inter-comparisons  with 
statistics  obtained  from  the  United  Kingdom  Meteorological  Office  (UKMO)  unified  model;  (c) 
investigate  the  scalability  of  COAMPS  on  ICM  computers,  a  200-processor  OPTERON  cluster  and  a 
Cray  XI;  and  (d)  develop  a  data  assimilation  scheme  that  can  assimilate  remotely-sensed  and  non- 
conventional  data  sources  with  a  special  emphasis  on  Doppler  radar  data.  Meeting  these  objectives 
will  allow  the  Polish  National  Air  Defense  to  issue  1-5  day  mesoscale  weather  forecasts  in  the  regions 
of  their  interest,  including  Poland  and  Central  Europe. 

APPROACH 

Our  approach  is  to  utilize  NOGAPS  for  initial  and  lateral  boundary  conditions,  and  COAMPS  for 
mesoscale  atmospheric  forecasts.  The  NOGAPS  fields  are  obtained  from  the  Global  Ocean  Data 
Assimilation  Experiment  (GODAE)  server  at  the  Elect  Numerical  Meteorology  and  Oceanography 
Center  (FNMOC)  and  transferred  to  ICM  in  automated  (machine-controlled)  efficient  and  stable  way, 
thanks  to  support  from  FNMOC.  We  ported  the  COAMPS  system  to  the  ICM  Einux  cluster  and  Cray 
XI  computers  and  measured  system  performance  and  scalability  using  tools  developed  in-house  for 
model  verification.  The  unique  aspect  of  our  capabilities  is  that  we  concurrently  run  the  UKMO 
mesoscale  model  on  a  grid  that  is  similar  to  the  one  used  by  COAMPS.  We  will  also  investigate  the 
time  evolution  of  the  conditional  forecast  (background)  error  probability  density  function  using  an 
ensemble  of  the  model  forecast  to  generate  background  error  statistics.  This  helps  us  to  identify  and 
understand  the  uncertainties  in  high-resolution  NWP  forecasts  on  high-impact  weather,  particularly 
extreme  rainfall.  Finally,  we  will  study  observational  error  characteristics  of  radar  reflectivity  and 
radar  radial  winds.  Such  observations  have  the  potential  to  provide  detailed  information  to  improve 
mesoscale  analyses  and  forecasts.  We  will  study  historical  weather  events  for  which  we  have  radar 
data  to  understand  the  observational  error  characteristics.  We  will  investigate  how  they  may  be 
applied  in  data  analyses  used  for  assimilating  radar  data  into  numerical  weather  prediction  models. 

WORK  COMPLETED 

During  FY07  we  accomplished  the  following  tasks:  (a)  made  a  revised  set-up  of  the  US  Navy 
CO  AMPS  system  on  ICM  machines  at  Warsaw  University  for  the  purpose  of  providing  operational 
support  to  the  general  public  using  the  ICM  meteorological  web  page,  (b)  improved  some  parts  of  the 
COAMPS  system  (more  frequent  boundary  conditions,  variable  step  in  interpolated  boundary  data),  (c) 
improved  the  optimal  interpolation  analysis  of  sea  surface  temperature  and  sea  ice  data  by  introducing 
sample  dependent  correction  of  the  background  error  covariances,  (d)  investigated  the  correlation  of 
errors  in  radar  reflectivity  data  using  historical  data  from  Baltex  Radar  Data  Centre,  (e)  worked  on  the 
development  of  the  entity  based  approach  to  the  verification  of  precipitation  patterns  predicted  by  the 
COAMPS  system,  and  (f)  worked  on  the  development  of  the  ensemble  Kalman  filter  approach  for  the 
purpose  of  assimilating  radar  reflectivity  and  radial  wind  data  in  COAMPS. 

RESULTS 

During  FY07,  we  sufficiently  increased  the  frequency  of  lateral  boundary  conditions  transmitted  from 
the  GODAE  server  to  our  model  domains.  Using  the  COAMPS  restart  capability,  we  implemented 


2 


three  restart  runs  with  different  frequeneies  for  the  interpolation  of  lateral  boundary  eonditions  (LBC). 
For  the  foreeast  period  from  0  to  24h  we  used  LBC  data  in  3h  intervals,  for  the  period  from  24  to  72h 
we  used  LBC  data  in  6h  intervals,  and  for  the  period  from  72  to  120h  we  used  LBC  data  in  12h 
intervals.  These  ehanges  in  LBC  frequeney  resulted  in  improved  predietion  of  seleeted  meteorologieal 
elements,  sueh  as  more  aeeurate  predietion  of  the  position  and  intensity  of  low  and  high  pressure 
eenters,  and  more  aeeurate  predietion  of  preeipitation  patterns.  We  also  introdueed  some  ehanges  into 
the  eonfiguration  of  our  model  domains.  Two  independent  domains  of  our  foreeast  were  designed. 

The  first  one,  ealled  NAE,  eovers  the  North  Atlantie  and  European  area,  and  eonsists  of  3  nested  grids 
with  different  spatial  resolutions  and  foreeast  lead  times.  The  eoarse  grid  has  193x127  grid  points  with 
a  mesh  size  of  39km  and  lead  time  of  120h;  the  seeond  one,  plaeed  over  Central  Europe,  has  169x217 
grid  points  with  a  mesh  size  of  39  km  and  lead  time  of  72h;  and  the  third  one,  plaeed  over  Poland,  has 
193x175  grid  points  with  mesh  size  of  4.3  km  and  lead  time  of  48h.  Our  seeond  domain,  ealled  ME, 
eovers  Middle  Eastern  eountries,  and  eonsists  of  two  nested  grids.  The  eoarse  grid  has  117x71  grid 
points  and  a  mesh  size  of  45  km,  and  the  inner  grid  has  259x127  grid  points  and  a  mesh  size  of  15  km. 
Both  grids  have  lead  times  of  72h.  This  domain  is  used  to  support  military  aetivities  of  Poland  in 
NATO  and  UN  operations. 

The  offieial  web  page  of  the  COAMPS  model  at  ICM  has  been  rewritten  and  tailored  to  the  needs  of 
our  general  publie  users.  The  page  eonsists  of  two  types  of  presentations:  maps  of  seleeted 
meteorologieal  elements  over  the  seeond  (13  km)  grid  of  our  NAE  domain  presented  in  3h  intervals 
until  72h,  and  the  meteograms  for  seleeted  eities  within  the  same  region.  The  user  has  the  option  to 
ehoose  any  point  in  the  region.  An  example  of  sueh  a  meteogram  is  shown  in  Eig.  1 .  The  time  series 
in  Ih  intervals  shows  the  temperature,  humidity,  rain,  mean  sea  level  pressure,  wind  speed  and 
direetion,  eloud  vertieal  extent,  and  amount  of  eloudiness  in  three  levels  (high,  medium,  and  low 
elouds)  in  separate  panels  on  the  diagram.  This  presentation  of  the  results  of  the  numerieal  weather 
foreeasts  is  very  popular  among  our  users. 


3 


sunrise;  06:47  CEST 
sunset:  t8:38  CEST 


16‘’56’E,  52«22'N 
Cx=7S,  y=132) 


height  (min,  centre,  ma>:) 
81  -  86  -  96  m 


meteo@icnn.edu.pl 


(C)  2007  ICM,  Uniwersytet  Warszawski 


Fig.  1.  Example  of  the  meteogram  constructed  from  COAMPS  forecasts 
for  the  13  km,  Central  Europe  grid.  Data  are  presented  in  Ih  intervals. 


In  addition  to  this  official  web  page  (http://new.meteo.pl/index_eng.php),  we  have  also  developed  a 
second  web  page  (http://coamps.icm.edu.pl),  where  selected  results,  both  at  the  near  surface  level  and 
at  few  levels  in  free  atmosphere  are  presented  for  all  nested  grids  for  both  the  NAE  and  ME  domains. 

During  EY07,  we  worked  also  on  improvements  to  the  COAMPS  data  assimilation.  The  COAMPS 
analysis  is  based  on  the  multivariate  optimum  interpolation  (MVOI)  analysis  scheme  described  in 
Goerss  and  Phoebus  (1992)  and  Barker  (1992).  The  MVOI  technique  uses  observational  data  to 
compute  increments  for  the  first-guess  fields  from  either  COAMPS  or  NOGAPS.  The  analysis 
variables  for  the  MVOI  are  geopotential  height,  and  the  u  and  v  wind  components.  The  first  guess 
fields  are  adjusted  based  on  observational  data  via  a  MVOI  analysis.  The  analysis  computes  height 
and  momentum  increments  based  on  differences  between  the  first  guess  fields  and  the  observations. 
The  COAMPS  system  also  has  an  ocean  analysis  component  -  the  NRE  Coupled  Ocean  Data 
Assimilation  (NCODA)  system  (Cummings,  2003).  We  utilize  the  2-dimensional  option  in  NCODA 
to  provide  updated  sea  surface  temperature  and  sea  ice  concentration  for  the  COAMPS  lower  boundary 
conditions.  The  analysis  background  is  generated  from  the  previous  NCODA  analysis.  NCODA 
computes  corrections  to  the  first-guess  fields  using  all  the  observations  that  have  become  available 
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since  the  last  analysis  was  made.  A  eomplete  derivation  and  deseription  of  the  MVOI  method  used  in 
NCODA  is  provided  in  Daley  (1991).  Essentially,  the  MVOI  problem  is  formulated  in  NCODA  as 

X  =  X  +  P^H"(HP,H"+  R)  -  '[y  -  H(x^]  (1) 

where  is  the  analysis,  Xh  is  the  baekground,  Pb  is  the  baekground  error  eo-varianee,  H  is  forward 
operator,  R  is  observation  error  eo-varianee,  and  j  is  the  observation  veetor.  The  observation  veetor 
eontains  all  of  the  temperature  and  sea  iee  observations  that  are  within  the  geographie  and  time 
domains  of  the  foreeast  model  and  update  eyele. 

A  forward  model  is  a  method  of  eonverting  a  foreeast  grid  variable  to  an  observed  variable.  The 
forward  operator  in  NCODA  is  simply  a  spatial  interpolation  of  the  foreeast  model  grid  to  the 
observation  loeation  performed  in  three  dimensions.  Thus,  HP^H^  is  ^proximated  direetly  by  the 
baekground  error  eo-varianee  between  observation  loeations,  and  PbET  direetly  by  the  error  eo- 
varianee  between  observation  and  grid  loeation.  The  quantity  [y-/T(jCb)]  is  referred  as  the  innovation 
veetor,  |j-/T(jCa)]  is  the  residual  veetor,  and  x^  -  Xb  is  the  inerement  veetor. 

The  baekground  error  eo-varianees  are  separated  into  a  baekground  error  varianee  and  eorrelation.  The 
eorrelation  Cb  is  further  separated  into  a  horizontal  Ch  and  a  vertieal  Cv  component.  All  correlations 
are  modelled  as  seeond  order  autoregressive  (SOAR)  funetions. 

The  model  error  parameterisation  imply  a  parameterisation  of  the  innovation  eo-varianee.  We  assume, 
that  we  have  a  eo-varianee  model  for  d,  whieh  has  a  number  of  unknown  parameters 

E(ddVS(a)  (2) 

One  way  to  calibrate  the  eo-varianee  model  iS'(a)  to  the  innovation  sample  d  is  to  estimate  a  using  the 
maximum  likelihood  method.  Eor  this  purpose,  we  assume  that  d  ~  N[0,  S(a)]  for  some  a  =  a*  and 
this  means  that  the  innovation  is  normally  distributed,  zero-mean  random  vector  whose  covarianee  is 
the  matrix  S(a*).  The  eonditional  probability  density  p(</|  a)  of  this  random  veetor  is  given  by  the 
Gaussian  density  funetion.  Given  the  innovation  sample  d,  the  maximum  likelihood  parameter 
estimate  of  a*  is  obtained  by  finding  the  value  of  a  for  whieh  the  probability  density  p(</|  a)  attains  a 
maximum.  This  is  equivalent  to  the  minimisation  of  the  maximum  likelihood  funetion  obtained  from 
natural  logarithm  of  the  probability  density. 

Our  experiments  eoneentrated  on  the  estimation  of  the  sealar  parameter  a  interpreted  as  a  sealing  faetor 
whieh  replaee  the  fixed  eo-varianee  matrix  So  by  the  parameterised  eo-varianee  aSg.  We  improve  the 
eonsisteney  between  the  eovarianee  model  used  in  our  seheme  and  the  actual  data  by  multiplying  the 
matrix  So  by  the  scaling  parameter  a.  In  our  case,  the  analysis  is  made  in  selected  sub-domains 
(volumes),  so  we  implemented  the  method  to  eaeh  of  volumes  separately.  Eor  the  grid  over  North 
Atlantie,  we  had  850  sueh  volumes,  for  Central  Europe  seas  we  used  180  volumes.  The  scaling  factor 
can  be  used  also  as  a  diagnostie  tool,  sinee  it  must  be  elose  to  unity  for  the  eovarianee  model  to  be 
eorreet.  Examining  a  in  the  North  Atlantic  region,  we  found  serious  errors  in  data  (where  a=9.3980), 
loeally  destroying  results  of  the  analysis  in  the  vieinity  of  this  point.  In  eontrol  runs,  these  erroneous 
data  were  also  assimilated,  however,  as  the  eo-varianee  model  did  not  follow  the  data  elosely,  the 
impaet  of  this  error  on  the  SST  analysis  was  smaller. 
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Fig.  2.  Comparison  of  SST  errors  in  the  operational  (left  panel)  and 
modified  (right  panel)  version  of  data  analysis  for  grid  2,  Central  Europe  region. 


By  examining  the  estimated  SST  fields  for  calibrated  and  control  runs,  we  find  that  our  approach 
produced  smoother  and  more  realistic  analyses.  In  the  case  of  quality  controlled  initial  data,  the  errors 
of  the  analysis  are  smaller  than  those  used  operationally.  However,  the  method  is  sensitive  to  the  gross 
errors  included  in  data,  so  quality  control  issues  for  this  method  are  essential. 

Another  field  of  activity  within  this  project  in  FY07  was  the  development  of  a  entity  based  scheme  of 
verification  of  precipitation  patterns  produced  by  the  COAMPS  system.  From  Fig.  3,  it  is  evident  that 
precipitation  patterns  estimated  using  parameterized  convection  (left  panel)  and  explicit  moist  physics 
(right  panels)  differs  in  position  and  intensity.  To  examine  the  behavior  of  the  model,  we  decided  to 
work  on  the  Ebert-McBride  technique  (EMT)  for  quantitative  precipitation  verification.  The  aim  of 
the  EMT  approach  is  to  verify  to  what  extent  the  forecast  entity  has  the  same  location,  shape,  and 
magnitude  as  the  observed  one,  with  resulting  error  statistics  based  on  the  properties  of  each  entity. 
This  type  of  verification  utilizes  a  pattern-matching  technique  to  horizontally  translate  the  forecast 
entity  over  the  observed  one.  The  best  match  can  be  determined  in  a  number  of  ways,  usually  by 
maximizing  the  correlation  coefficient  or  by  minimizing  the  total  error. 
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Fig.  3.  Precipitation  forecasts  at  13  km  grid  (left  panel) 
and  4  km  grids  (right  panels)  in  moving  frontal  situation. 


We  have  now  developed  software  whieh  enables  us  to  produee  the  forecasts  and  radar  observation  on 
common  projections  and  with  the  same  resolutions.  The  verification  area  is  limited  to  the  area  covered 
by  the  radars  (in  our  case,  the  Baltic  Sea  catchment).  We  also  are  able  to  divide  the  forecasts  and 
observations  into  separate  entities  and  our  algorithm  keeps  track  of  the  evolution  of  these  entities.  We 
are  working  on  statistics  such  as  the  mean  horizontal  displacement  of  the  forecast,  the  error  in  the 
forecast  and  observed  rain  area,  main  and  maximum  rain  rates,  rain  volume,  and  the  pattern  correlation 
of  the  corrected  forecast.  Although  the  position  of  the  forecast  patterns  is  still  different  from  the 
observed  one,  the  main  features  of  the  convective  systems  are  simulated  reasonably  well.  The  example 
presented  on  Fig.  4  comes  from  our  4  km  run,  in  which  explicit  moist  physics  is  used.  The  cross- 
section  shows  well  developed  cumulonimbus  clouds  with  an  anvil  top. 

Another  field  of  activity  within  this  project  in  FY07  was  the  work  on  an  ensemble  data  assimilation 
system  based  on  the  Kalman  square  root  filters.  Ensemble  data  assimilation  systems  transform  a 
forecast  ensemble  into  an  analysis  ensemble  with  appropriate  statistics.  This  can  be  done  statistically 
by  treating  the  observations  as  random  variables,  or  deterministically,  requiring  that  the  covariance  of 
the  updated  ensemble  satisfy  the  Kalman  filter  analysis  error  covariance  equation.  We  implemented 
the  deterministic  analysis  ensemble  update,  which  is  a  Monte  Carlo  implementation  of  the  Kalman 
square  root  filters.  Ensemble  square  root  filters  are  not  unique,  since  different  ensembles  can  have  the 
same  covariance.  This  non-uniqueness  had  led  to  the  development  of  several  different  algorithms  for 
updating  the  analysis  ensemble.  In  our  implementation  of  the  square  root  filter,  we  used  the  Whitaker 
and  Hamil  (2002)  method. 
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Fig.  4.  Vertical  cross-section  trough  the  frontal  system.  4.3  km  grid  over  Poland,  22  Jul  2007. 


An  additional  processing  of  the  ensemble  covariances  is  introduced  to  avoid  the  filter  divergence 
caused  by  sampling  errors.  The  techniques  commonly  used  (and  implemented  in  our  work)  are 
distance  dependent  covariance  localization  and  covariance  inflation.  Covariance  localization  is  a  filter 
that  forces  the  ensemble  covariances  to  go  to  zero  at  some  horizontal  distance  L  from  the  observation 
being  assimilated.  Covariance  inflation  simply  inflates  the  deviations  from  the  ensemble  mean  first 
guess  by  some  factor  greater  than  1 .0  for  each  member  of  the  ensemble,  before  the  computation  of  the 
background-error  covariances  and  before  any  observations  are  assimilated. 

PERSONNEL  EXCHANGES  AND  TRAVEL  COMPLEMENTED 

Marcin  Witek,  University  of  Warsaw  -  visited  NRL,  Monterey,  CA  to  work  with  Piotr  Flatau  on 
development  of  aerosol  parameterization  in  NAAPS  and  COAMPS  models. 

Richard  Hodur,  NRL,  Monterey,  CA  -  in  November  2006  visited  Warsaw  University  for  two 
lectures  on  recent  developments  in  numerical  weather  prediction  at  NRL  (one  lecture  was  presented  at 
ICM  and  the  second  at  Military  Meteorological  Office)  and  working  discussions  on  COAMPS  issues. 
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Bogumil  Jakubiak,  University  of  Warsaw  -  participated  in  EUMET/SRNWP  Workshop  on  “High 
resolution  data  assimilation  with  emphasis  on  the  use  of  moisture-related  observations”.  21-23  Mareh 
2007,  Norrkoping,  Sweden. 

Oskar  Kapala,  University  of  Warsaw  -  in  April  2007  visited  NRL,  Monterey,  CA,  to  work  with 
Keith  Sashegyi  on  nested  3D-VAR  algorithm  and  to  discuss  with  Rich  Hodur,  Piotr  Elatau  and  Maria 
Elatau  the  new  developments  in  model  dynamies.  He  discussed  also  with  Ted  Holt  the  possibility  to 
implement  a  new  land-surfaee  seheme  into  COAMPS  model. 

Bogumil  Jakubiak,  University  of  Warsaw  -  participated  in  EGU  General  Assembly,  Vienna  Austria 
15-20  April  2007  giving  one  oral  and  two  poster  presentations. 

IMPACT/APPLICATIONS 

The  operational  version  of  the  high  resolution  mesoseale  model  was  implemented.  This  will  improve 
5-day  forecasting  for  aviation  and  for  support  of  Polish  troops  in  Afghanistan  and  Iraq. 

TRANSITIONS 

None. 

RELATED  PROJECTS 

US  GODAE  project  -  The  Global  Oeean  Data  Assimilation  Experiment  (GODAE)  provides  regular, 
eomplete  deseriptions  of  the  state  of  the  oeean  and  the  atmosphere  in  support  of  operational 
oeeanography  and  oceanographie  researeh.  We  use  the  GODAE  server  data  for  our  COAMPS 
predietions.  The  GODAE  Monterey  server  is  maintained  by  ENMOC  and  is  sponsored  by  the  Office 
of  Naval  Researeh  (ONR). 

COST  Action  731  project  -  Propagation  of  uncertainty  in  advanced  meteo-hydrological  forecast 
systems.  Within  this  action,  we  started  to  develop  a  radar  data  assimilation  scheme  using  the  ensemble 
Kalman  filter  approach. 
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